---
title: "RCT"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

library(dplyr)
library(sandwich)
library(lmtest)
library(MASS)
library(emmeans)
library(readr)
library(ggplot2)
library(effects)
library(tidyr)
library(BalanceR)
library(memisc)

dat <- read_csv("data/Iraq_election2025.csv")

# RCT 
dat <- dat |> 
  mutate(
    Voting_pre  = as.numeric(Q8),
    Voting_post = as.numeric(Q11),
    stimulus = as.numeric(Stimulus),
    Stimulus = factor(stimulus, level = c(1, 2),
                      labels = c("Control", "Treatment"))
  ) |> 
  filter(!is.na(Voting_pre), !is.na(Voting_post), !is.na(Stimulus)) |> 
  mutate(Gender = D1,
      Age = D2 + 17,
      Education = D3,
      Income = D5,
      Sect = case_when(D6 == 1 ~ "Sunni",
                       D6 == 2 ~ "Shia",
                       D6 == 3 ~ "Christian",
                       D6 == 4 ~ "Kurd",
                       D6 == 5 ~ "Other"),
      Shia = if_else(D6 == 2, 1, 0),
      Support_same_ethno_sectarian_origin = Q14_3)

# delta
dat <- dat |> 
  mutate(delta = Voting_post - Voting_pre)

```

# balance check
```{r}
table(dat$Stimulus, dat$Voting_pre)
table(dat$Stimulus, dat$Voting_post)

chisq.test(table(dat$Stimulus, dat$Voting_pre))

dat <- dat |> 
  mutate(sect = D5)
balance <- BalanceR(data = dat, 
                    group = Stimulus,
                    cov = c(Gender, Age, Education, Income, sect))
print(balance, digits = 3)
summary(balance, digits = 5)
plot(balance, point.size = 5, text.size = 18)

```

# Difference Model or Change Score Model, as rubustness check
```{r}
# model 9 and 10
reg20 <- lm(delta ~ Stimulus, data = dat)
coeftest(reg20, vcov = vcovHC(reg20, type = "HC1"))
summary(reg20)

reg21 <- lm(delta ~ Stimulus + Gender + Age + Education + Income, data = dat)
summary(reg21)

# Figure 1: plot the difference of average
delta_summary <- dat |> 
  group_by(Stimulus) |> 
  summarise(
    mean_delta = mean(delta, na.rm = TRUE),
    sd_delta = sd(delta, na.rm = TRUE),
    n = n(),
    se_delta = sd_delta / sqrt(n)
  ) |> 
  mutate(
    lower = mean_delta - 1.96 * se_delta,
    upper = mean_delta + 1.96 * se_delta
  )

ggplot(delta_summary, aes(x = Stimulus, y = mean_delta, fill = Stimulus)) +
  geom_bar(stat = "identity", width = 0.6, color = "black") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  ylab("Mean Change in Voting Intention (Δ)") +
  xlab("Stimulus") +
  theme_minimal() +
  scale_fill_manual(values = c("lightblue", "salmon"))
```

# OLS: ANCOVA
```{r}
# Model 1
reg01 <- lm(Voting_post ~ Stimulus + Voting_pre, data = dat)
coeftest(reg01, vcov. = vcovHC(reg01, type = "HC1"))
summary(reg01)

reg02 <- lm(Voting_post ~ Stimulus * Voting_pre, data = dat)
summary(reg02)

# Model 2
reg03 <- lm(Voting_post ~ Stimulus + Voting_pre +
              Gender + Age + Education + Income, data = dat)
summary(reg03)

# Model 3
reg04 <- lm(Voting_post ~ Stimulus * Support_same_ethno_sectarian_origin + Voting_pre, data = dat)
summary(reg04)

reg05 <- lm(Voting_post ~ Stimulus * Shia + Voting_pre, data = dat)
summary(reg05)

# Model 4
reg06 <- lm(Voting_post ~ Stimulus * Support_same_ethno_sectarian_origin + Voting_pre +
              Gender + Age + Education + Income, data = dat)
summary(reg06)

# Figure 4: marginal effect of reg04
eff <- effect("Stimulus:Support_same_ethno_sectarian_origin", reg04)
eff_df <- as.data.frame(eff)

ggplot(eff_df, aes(x = Support_same_ethno_sectarian_origin, 
                   y = fit, color = Stimulus)) +
  geom_line(size = 1.2) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = Stimulus), alpha = 0.15, color = NA) +
  labs(
    x = "Degree of importance placed on leaders of the same sectarian/ethnic origin in party choice",
    y = "Predicted voting intent（post stimulus）",
    title = "Interaction effect: Importance of same origin × Treatment (Stimulus）"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"))
```

# regression table
```{r}
table <- mtable("Model 1" = reg01, "Model 2" = reg03,
                "Model 3" = reg04, "Model 4" = reg06,
                "Model 5" = reg20, "Model 6" = reg21,
               summary.stats=c("adj. R-squared","F","p", "N"))
print(table)
write_html(table, file = "result/RCT_regression.html")

table1 <- mtable("Model 1" = reg01, "Model 2" = reg03,
                 "Model 3" = reg04, "Model 4" = reg06,
                summary.stats=c("adj. R-squared","F","p", "N"))
print(table1)
write_html(table1, file = "result/Table1.html")

app_table1 <- mtable("Model 5" = reg20, "Model 6" = reg21,
                    summary.stats=c("adj. R-squared","F","p", "N"))
print(app_table1)
write_html(app_table1, file = "result/App_table1.html")

```

# Ordinal as robustness in appendix
```{r}
polr1 <- MASS::polr(factor(Voting_post) ~ Stimulus + factor(Voting_pre),
                    data = dat, Hess = TRUE)
summary(polr1)

polr2 <- MASS::polr(factor(Voting_post) ~ Stimulus * Support_same_ethno_sectarian_origin 
                    + factor(Voting_pre) + Gender + Age + Education + Income,
                    data = dat, Hess = TRUE)
summary(polr2)


library(broom)
library(modelsummary)

tidy_polr1 <- tidy(polr1)
tidy_polr2 <- tidy(polr2)

models_polr <- list(
  "Ordered Logit" = polr1,
  "Ordered Logit + Interaction" = polr2
)

modelsummary(
  models_polr,
  statistic = "({std.error})",
  stars = TRUE,
  coef_map = c(
    "StimulusTreatment" = "Treatment",
    "Support_same_ethno_sectarian_origin" = "Ethno-sectarian importance",
    "StimulusTreatment:Support_same_ethno_sectarian_origin" =
      "Treatment × Ethno-sectarian importance",
    "1|2" = "Cutpoint 1|2",
    "2|3" = "Cutpoint 2|3",
    "3|4" = "Cutpoint 3|4"
  ),
  gof_map = data.frame(
    raw = c("nobs", "AIC"),
    clean = c("N", "AIC"),
    fmt = c(0, 2)
  ),
  output = "result/appendix_ordered_logit_full.html"
)

```

# statistics table
```{r}
dat_desc <- dat |>
  dplyr::select(
    Voting_pre,
    Voting_post,
    stimulus,
    delta,
    Support_same_ethno_sectarian_origin
  )

psych::describe(dat_desc)[, c("n", "mean", "sd", "min", "max")]
```

